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Abstract 

In this Letter, we present results of an extensive Monte Carlo study of the O/W(110) 
system under non-equilibrium conditions. We study the mean square displacements 
and long wavelength density fluctuations of adatoms. From these quantities, we 
define effective and time-dependent values for the collective and tracer diffusion 
mobilities. These mobilities reduce to the usual diffusion constants when equilibrium 
is reached. We discuss our results in view of existing experimental measurements 
of effective diffusion barriers, and the difficulties associated with interpreting non- 
equilibrium data. 
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1 Introduction 



Surface diffusion plays a significant role in many surface phenomena such as 
epitaxial growth, catalysis, and ordering [1]. These phenomena are interest- 
ing both from the point of view of basic research and applied science. Many 
new experimental techniques [2-4] have been developed for the measurement 
of surface diffusion either under equilibrium conditions or slight deviations 
from equilibrium. There is a corresponding increase in our theoretical under- 
standing of this important quantity. Under near equilibrium conditions, the 
theoretical description of surface diffusion is well established based on the lin- 
ear response theory by Kubo [5]. Equivalently, the problem can be described 
by the diffusion equation governing the evolution of the density profile [6]. 
However, in many non-equilibrium situations such as surface growth, there is 
no unique way of defining a diffusion constant D. Yet the diffusive motion of 
the adatoms clearly plays an important role in such processes. In this work, we 
will introduce the concepts of tracer and collective mobilities which character- 
ize the motion of adatoms under non-equilibrium situations. These quantities 
are defined as generalizations of the usual diffusion constants to which they 
reduce in the appropriate limits. To study the non-equilibrium mobilities, we 
have undertaken an extensive Monte Carlo study of the O/W(110) system 
using a lattice-gas model. We consider the non-equilibrium ordering dynamics 
of this system after a quench from a totally random initial state to a tempera- 
ture characterizing an ordered state. We then divide the time scales into slices 
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according to the decay of the excess energy. In each time slice, we introduce 
a definition for tracer and collective mobility. These mobilities are then fitted 
to an Arrhenius form to extract the effective diffusion barriers that are shown 
to be strongly time- dependent during the ordering process. For comparison, 
the equilibrium properties have also been determined. We discuss our results 
in light of the experimental data by Tringides et al. [7-10], and consider the 
difficulties associated with interpreting non-equilibrium measurements. 

2 Model and methods 

The phase diagram of the O/W(110) system is fairly well determined through 
a series of experimental studies [9,11,12]. Its main features have also been 
obtained from theoretical calculations with a lattice-gas Hamiltonian including 
pair interactions up to fifth neighbors and triplet interactions [13]. We use this 
lattice-gas Hamiltonian for our present study of non-equilibrium properties of 
this system through Monte Carlo (MC) simulations. For this work, we focus 
at the coverage 9 = 0.45. At this coverage, the experimental value for the 
critical temperature of the order-disorder transition is T c cxp ps 710 K [12]. 
Above T c , the adsorbate is in a disordered phase. Below T c , the system is 
characterized by an ordered p(2 x 1) or equivalently the p(l x 2) phase. In 
this work, we study the ordering dynamics of the system by starting from an 
initial state, which is completely random such as that obtained from random 
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deposition at low temperatures. The system is then allowed to evolve towards 
the equilibrium ordered phase at various temperatures ranging from 0.655 T c 
to 0.893 T c . This procedure corresponds to instantaneous heating in an actual 
experimental situation. 

To perform our MC simulations, we chose the transition dynamics algorithm 
(TDA) [14]. Within TDA the transition rate Wij, from an initial state % with 
energy to a final state / with energy Ef, is decomposed into two steps by 
introducing an intermediate state I with energy Ej = (Ei + Ef)/2 + A, where 
the quantity A characterizes the activation barrier in the zero coverage limit 
due to the substrate-adatom interaction. The rate w it f is then the product of 
the two rates Wij and wij, which are taken to be of the standard Metropolis 
form [15]. The TDA describes the classical diffusion barrier more realistically 
than other transition rate algorithms [14,15], in which the effect of the saddle 
point of the adiabatic surface potential is not taken into account. 

For the present study, we chose A = 0.0437 eV. This value is believed to be 
much lower than the true value which should be closer to the experimentally 
observed barrier of 0.5 to 0.6 eV [9,16] in the disordered phase. Our choice is 
necessitated by the need to speed up the jump rate in the numerical simula- 
tions at low temperatures. We have done some tests and found that the effect 
of A and the adatom-adatom contributions to the diffusion barrier are approx- 
imately additive. Therefore the barriers calculated here should be increased by 
about 0.5 to 0.6 eV for comparison with the experimentally observed values. 
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To analyse the ordering dynamics in terms of quasi-equilibrium concepts, we 
have to divide the total ordering period into different time regimes. These slices 
of time obviously have to be different at different temperatures because of the 
change in the rate of ordering. For this purpose, we first calculate the time- 
dependence of the excess energy of the system, namely E x (T,t) = E(T,t) — 
E(T, oo), after an instantaneous quench at time t — from a completely 
random state to a temperature T < T c . Here E(T,t) is the energy of the 
system at temperature T and time t. We then introduce the normalized excess 
energy F(T, t) = E X (T, t)/E x (T, 0), which has the maximum F(T, 0) = f. 
The equivalent time regimes at different temperatures are chosen as intervals 
between times t n (T) which satisfy F(T,t n {T)) = exp(— n), where the integer 
n > 0. Typically, we have used five such time slices here. 

Now we come to the definition of the non-equilibrium mobilities. In the case 
of tracer mobility, we consider the quantity 

1 N 

&t) = W J2(\R k a (tn(T)+5t)-R k a (tn(T))\ 2 ), (f) 

where a = x,y and the sum is over iV particles to improve statistics. The 
two independent spatial components (x, y) of the position vector R k (t) for a 
particle k at time t are denoted by R k (t), i.e. R k (t) = (R k (t),Ry(t)), and ( ) 
is used to denote configuration averaging. Within each slice of time, the time 
difference St obeys < St < t n+1 (T) — t n (T). The tracer mobility D^ aa is 
then defined as the effective slope of £^ within the given time regime n. In 
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equilibrium and for St — > oo, Eq. (1) is just the tracer diffusion constant D T 
times St. If the system approaches equilibrium as a function of time, ^ n ^/St 
tends towards D T) too. 

To study collective mobility, we introduce the time- dependent density-fluctuation 
autocorrelation function 

S^(R,R',St) = (5p(R,t n (T) + St)Sp(R',t n (T))), (2) 

— * — » 

where R and R denote position vectors on a lattice, and the density ffuctu- 

— * — * — * 

ations are given by Sp(R,t) = p(R,t) — (p(R,t)). Again, < St < t n+ i(T) — 
t n (T) within each time regime n. In the limit of long times (large n and St) 
and in the hydrodynamic regime, i.e. in equilibrium, the Fourier transform of 
this correlation function decays as S(k, St) = S(k, 0) exp(— k • D • kSt). In the 
non-equilibrium situation, we consider log S^- n \k, St) over the given regime n 
and define its effective slope divided by k 2 as the collective mobility D^ aa . 
In practice this method was carried out with separate sine and cosine trans- 
forms of the density fluctuations, following the approach in Ref. [17]. A more 
detailed presentation of the various methods used will be given elsewhere [18]. 

The MC calculations were carried out in a M x M lattice, the system size 
being M = 30. Although it is rather small, it is large enough when one is 
not close to the critical region. The number of independent samples varied 
between 1000 and 10000. 
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3 Results and Discussion 

We want to examine first the validity of extracting the mobilities from the 
quantities ^(St) and log S( n \k, St) by assuming their linear dependence on 
St within each time slice. A typical example of the behavior of £%£(5t) is given 
in Fig. 1. The deviations from linear behavior are not very large at any time, 
therefore our definition for is well justified. A similar result is also found 
for log S( n \k, St) used to extract . Furthermore, it is evident from Figs. 2 
(a) and 2 (b) that within each time slice n the temperature dependence of the 
resulting tracer and collective mobilities are well described by an Arrhenius 
form ~ exp(— /3Ea), with Ea denoting the effective activation barrier 
for diffusion. This allows us to determine the time-dependence of Ea during 
the ordering process. The results for Ea based on the tracer and collective 
mobilities are shown in Fig. 3. Three interesting features emerge. At early 
times, the adsorbate is in a disordered configuration and the adatom-adatom 
interaction contributions to the activation barrier largely cancel out. Thus, 
Ea is approximately just the intrinsic barrier A. At intermediate times, Ea 
increases rapidly, with the barrier for the tracer mobility always larger than 
the corresponding one for the collective mobility. Finally, at long times Ea 
approaches the equilibrium value E A q = 0.297 ±0.008 eV given by the dashed 
line. We note that, in the equilibrium, £U's for tracer and collective diffusion 
are indeed equal within error limits. 
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We now discuss the relevance of our results to some experiments on the 
O/W(110) system by Tringides et al. [7-10]. These authors studied the or- 
dering dynamics of the p(2 x 1) phase at 6 ~ 0.5 following an up-quench in 
temperature as in our simulation work. They found that the time-dependent 
average domain size L(t) followed the growth law L(t) = A(T)t^ with ps 0.28, 
where is the kinetic growth exponent. Using dimensionality arguments, they 
argued that the prefactor A(T) can be related to an effective diffusion con- 
stant D through the relation A(T) oc D^. Thus, the measured temperature 
dependence of A(T) allows the determination of an effective diffusion barrier 
E A . For this they obtained [7,9] the value of E A P = 0.61 ± 0.11 eV in con- 
trast to the equilibrium measurements of Gomer et al. [16,19] who obtained 
£,ex P ,eq _ ^ g _|_ g e y a t $ ~ o.56 in the ordered phase. We performed a 
similar analysis for the ordering process in our simulations. We used E~ l (T, t) 
as a measure for L(t) [20] and found that the growth law is valid over a certain 
period of time only; namely from the end of an initial transient period to the 
point when the finite size or pinning effects come into play. In our case, we 
obtain <fi m 0.5 in the regime corresponding to the time slices n = 2,3. This 
is illustrated in Fig. 4. We note that at large times, the growth exponent is 
not well defined due to finite-size effects that are rather pronounced after the 
regime n — 3. The value of the effective barrier Ea extracted from the relation 
A(T) oc D* is 0.192 ± 0.021 eV. Our time-dependent results for E A based on 
the mobilities for these intermediate time regimes range between 0.88 E A 
and 1.1 E A - Thus, we conclude that the procedure of using the dimensional 
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analysis to extract an effective barrier from the prefactor in the power law 
growth is indeed valid, and it represents the barrier for the dominant diffusion 
process in that time regime. Concerning the absolute values of the diffusion 
barriers, we have to add about 0.6 eV to them as discussed earlier to account 
for our choice of A. This brings our results for the effective barrier in the 
power law growth regime to 0.8 eV and the final equilibrium barrier to 0.9 eV. 
The equilibrium value is in agreement with the experimental value. Also we 
see the same that the adatom-adatom contribution is most dominant in the 
equilibrium ordered phase and the effective barrier measured during the or- 
dering processes has a smaller value than the equilibrium barrier [21,22]. The 
discrepancies between our calculated value for Ea and the measured value are 
not too surprising, since ~ 0.28 is considerably lower than the value ~ 0.5 
observed in our simulations. For the one-component p(2 x 1) phase whose or- 
der parameter is not conserved, the theoretically expected value for is 1/2 
[20]. Experimentally this value has been observed in the 0/W(112) system 
[23], for example. The experimental data for the O/W(110) system in Refs. 
[7,9] probably belongs to the early time rather than the intermediate time 
regime, since the maximum mean diameter of the growing domains during 
the ordering experiments is smaller than six lattice spacings [7,9]. The E^ p 
they measure though represents the effective diffusion barrier for a different 
configuration than the one in the intermediate time regime. The exact cause 
for the failure of the experiment to reach the t 1 / 2 growth regime is still unclear 
[10]. 
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We believe that in other non-equilibrium methods such as profile evolution 
techniques, useful information for an effective mobility is also best extracted 
from an intermediate time regime. Otherwise, the initial behavior with large 
deviations from equilibrium behavior would yield results difficult to interpret. 
Unlike the theoretical simulation studies where we can define the intermediate 
regimes rather precisely, the practical difficulties in actual non-equilibrium 
measurements would be to identify the proper intermediate time regimes. In 
the domain ordering dynamics, we have seen that the power law growth regime 
can be identified for this purpose. For other non-equilibrium situations, it is 
not clear how to establish experimentally similar criteria. 
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Fig. 1. Results for ^ (St) with five values of n, after a quench at t = from a 
completely random state to T = 0.714 T c . Time difference 5t is given in Monte 
Carlo steps (mcs). In the inner figure, a part of the small time behavior has been 
magnified. The results for £yy (not shown here) are similar. 



Fig. 2. (a) Results for the tracer mobilities Dj, xx as an Arrhenius plot. In addition 
to the equilibrium case that is given for the purpose of comparison (stars and dotted 
line), results up to n = 4 are presented. The results for Dj? yy (not shown here) are 
similar, (b) A similar plot in the case of collective mobility D^ xx . Due to separate 
sine and cosine transforms of the density fluctuations, two values for each pair of T 
and n are given. The results for (not shown here) are alike. 
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Fig. 3. Results for the effective activation energy Ea versus time. Results with five 
values of n based on the tracer and collective mobility are shown with full squares 
and open circles, respectively. Error bars are smaller then the size of the symbol for 
tracer and approximately of the same size for collective mobility. The equilibrium 
limit Ea/&a = 1 with E A = °- 297 eV is S iven W the dashed line. The fit for E A 
based on the A(T) data is given in the inset. (See text for details.) 



Fig. 4. A typical example of the time evolution of E~ data at T = 0.655 T c as a 
log-log plot. Time is given in Monte Carlo steps (mcs), and the system size used was 
120 x 120. A reference curve with = 1/2 (dotted line) over the regimes n = 2,3 
is also given. The kinetic growth exponent 4> and the A(T) data were evaluated 
between 6.9 and 8.3 for logt, yielding 0.514 for (f>. 
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Figure 1 by Vattulainen et al. 
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